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Abstract 

We numerically study dynamical behaviors of the quasiperiodically forced Hodgkin-Huxley neu- 
ron and compare the dynamical responses with those for the case of periodic stimulus. In the 
periodically forced case, a transition from a periodic to a chaotic oscillation was found to occur via 
period doublings in previous numerical and experimental works. We investigate the effect of the 
quasiperiodic forcing on this period-doubling route to chaotic oscillation. In contrast to the case 
of periodic forcing, new type of strange nonchaotic (SN) oscillating states (that are geometrically 
strange but have no positive Lyapunov exponents) are found to exist between the regular and 
chaotic oscillating states as intermediate ones. Their strange fractal geometry leads to aperiodic 
"complex" spikings. Various dynamical routes to SN oscillations are identified, as in the quasiperi- 
odically forced logistic map. These SN spikings are expected to be observed in experiments of the 
quasiperiodically forced squid giant axon. 
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I. INTRODUCTION 



To probe dynamical properties of a system, one often applies an external stimulus to the 
system and investigate its response. Particularly, periodically stimulated biological oscilla- 
tors have attracted much attention in various systems such as the embryonic chick heart-cell 
aggregates [1] and the squid giant axon O |3]. These periodically forced systems have been 
found to exhibit rich regular and chaotic behaviors P|. Recently, similar lockings and chaotic 
responses have also been found in neocortical networks of pyramidal neurons under periodic 
synaptic input [5J. On the contrary, quasiperiodically forced case has received little atten- 
tion [B]. Hence, intensive investigation of quasiperiodically forced biological oscillators is 
necessary for understanding their dynamical responses under the quasiperiodic stimulus. 

Strange nonchaotic (SN) attractors typically appear between the regular and chaotic 
attractors in quasiperiodically forced dynamical systems [71-[T8]. They exhibit some prop- 
erties of regular as well as chaotic attractors. Like regular attractors, their dynamics is 
nonchaotic in the sense that they do not have a positive Lyapunov exponent; like usual 
chaotic attractors, they have a geometrically strange fractal structure. Here, we are inter- 
ested in dynamical responses of neural oscillators subject to quasiperiodic stimulation. SN 
oscillations are expected to occur in quasiperiodically forced neural oscillators. 

This paper is organized as follows. In Sec. [TT| we study dynamical responses of the 
quasiperiodically forced Hodgkin-Huxley (HH) neuron which was originally introduced to 
describe the behavior of the squid giant axon [T^. As a dc stimulus passes a threshold 
value, a self-sustained oscillation (corresponding to a spiking state) is induced. Effect of 
periodic forcing on this HH oscillator was previously studied [3l EQ]. Thus, regular (such as 
phase locking and quasiperiodicity) and chaotic responses were found. We note that similar 
dynamical responses were also observed in experimental works of the periodically forced 
squid giant axon [21 E]- In this way, the model study on the HH neuron may be examined in 
real experiments of the squid giant axon. Here, we numerically study the case that the HH 
oscillator is quasiperiodically forced at two incommensurate frequencies and compare the 
dynamical responses with those for the periodically forced case. In the periodically forced 
case {i.e., in the presence of only one ac stimulus source), a transition from a periodic to 
a chaotic oscillation was found to occur via period-doubling bifurcations [3l |2^. Effect of 
the quasiperiodic forcing on this period-doubling route to chaotic oscillation is particularly 
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investigated by adding another independent ac stimulus source. Thus, unhke the case of 
periodic forcing, new type of SN oscillating states are found to occur between the regular 
and chaotic oscillating states as intermediate ones. As a result of their strange geometry, SN 
oscillating states give rise to the appearance of aperiodic complex spikings, as in the case of 
chaotic oscillations. Diverse dynamical routes to SN oscillations are identified, like the case 
of the quasiperiodically forced logistic map 0. It is also expected that these SN spikings 
might be observed in experimental works on the quasiperiodically forced squid giant axon. 
Finally, a summary is given in Sec. |III[ 



II. SN OSCILLATIONS IN THE QUASIPERIODICALLY FORCED HH OSCIL- 
LATOR 



We consider the conductance-based HH neuron model which serves as a canonical model 
for tonically spiking neurons. The dynamics of the HH neuron, which is quasiperiodically 
forced at two incommensurate frequencies /i and /2, is governed by the following set of 
differential equations: 

dV 

= -QNam^'hiV - V^a) " 9Kn\V - Vk) -gdV- Vl) + /e.t, (la) 

dx 

— = a^{y){l - x) - I3^{y)x] x = m,h,n, (lb) 

where the external stimulus current density (measured in units of /lA/cm^) is given by 
hxt = Idc + ^1 sin(27r/it) + A2 sin(27r/2t), Idc is a dc stimulus, Ai and A2 are amplitudes of 
quasiperiodic forcing, and ijj{= /2//1) is irrational (/i and /2: measured in units of kHz). 
Here, the state of the HH neuron at a time t (measured in units of ms) is characterized 
by four variables: the membrane potential V (measured in units of mV), the activation 
(inactivation) gate variable m (h) of the Na^ channel [i.e., the fraction of sodium channels 
with open activation (inactivation) gates], and the activation gate variable n of the 



channel {i.e., the fraction of potassium channels with open activation gates). In Eq. (la), C 
represents the membrane capacitance per surface unit (measured in units of /iF/cm^) and 
the total ionic current lion consists of the sodium current In a, the potassium current Ik, and 
the leakage current II- Each ionic current obeys Ohm's law. The constants g^a, Ok, and 
gi are the maximum conductances for the ion and the leakage channels, and the constants 



VNa, Vk, and Vl are the reversal potentials at which each current is balanced by the ionic 
concentration difference across the membrane. The three gate variables obey the first-order 



kinetics of Eq. (lb). The rate constants are given by 



aUV) = 0.07eM-(V - K)/20|. ^V) = ,,^1(30 _ (^^ ;.^)}/io] + l ' 

""O') = exp[{To-(^7im-i - = - P=) 

where Vr is the resting potential when lext = 0. For the squid giant axon, typical values 
of the parameters (at 6.3 °C) are [21]: C = l/xF/cm^, g^a = 120mS/cm^, = 36mS/cm^, 
gL = 0.3mS/cm2, Vjva = 50 mV, Vk = -77 mV, Vl = -54.4 mV, and K = -65 mV. 

To obtain the Poincare map of Eq. ([T]), we make a normalization fit — )■ t, and then 
Eq. ([T]) can be reduced to the following differential equations: 
dV 



dt 



Fi(x,( 
1 , 



-gNani^h{V - Vno) - 9Kn\V - Vk) - g^V - Vl) + hxt], (3a) 



C fi 

(iTfl 1 

— = F2{^,9) = -[am{V){l-m)-P^{V)ml (3b) 
^ = F,{^,d) = j[an{V){l-h)-(5^{V)hl (3c) 

(ITI 1 

- = F4(x,0) = -K(\/)(l-n)-/3„(\/)n], (3d) 



de 

dt 



u (mod 1), (3e) 



where x[= (xi, X2, Xs, X4)] = {V,m,h,n) and /ext = -^dc + ^1 sin(27rt) + sin(27r6'). The 
phase space of the quasiperiodically forced HH oscillator is six dimensional with coordinates 
V, m, h, n, 9, and t. Since the system is periodic in 9 and t, they are circular coordinates 
in the phase space. Then, we consider the surface of section, the V-m-h-n-9 hypersurface 
a.t t = n {n: integer). The phase-space trajectory intersects the surface of section in a 
sequence of points. This sequence of points corresponds to a mapping on the five-dimensional 
hypersurface. The map can be computed by stroboscopically sampling the orbit points v„ 
[= (x„, 9n)] at the discrete time n (corresponding to multiples of the first external driving 
period Ti). We call the transformation v„ — )■ v,„+i the Poincare map, and write v^+i = 

Piyn). 



Numerical integration of Eqs. ([T| and ^ is done using the fourth-order Runge-Kutta 
method. Dynamical analysis is performed in both the continuous-time system {i. e., flow) and 
the discrete-time system {i.e., Poincare map). For example, the time series of the membrane 
potential V{t) and the phase flow are obtained in the flow. On the other hand, the Lyapunov 
exponent [22j and the phase sensitivity exponent |12j of an attractor are calculated in the 
Poincare map. To obtain the Lyapunov exponent of an attractor in the Poincare map, 
we choose 20 random initial points {(1^(0), mj(0), /;,j(0), ^^(O), 6'j(0)); z = 1,...,20} with 
uniform probability in the range of V,{0) G (-60,0), mi{0) G (0.1,0.9), hi{0) G (0.1,0.2), 
nj(0) G (0.5,0.7) and 6'j(0) G [0,1). For each initial point, we get the Lyapunov exponent, 
and choose the average value of the 20 Lyapunov exponents. (The method of obtaining the 
phase sensitivity exponent will be explained below.) 

In the presence of only the dc stimulus {i.e., A\ = A2 = 0), a transition from a resting 
state to a periodic spiking state occurs for J^c = /^c (— ^-^"^ /iA/cm^) via a subcritical 
Hopf bifurcation when the resting state absorbs the unstable limit cycle born via a fold 
limit cycle bifurcation for J^c — 6.26 /lA/cm^ [221 121] • Thus, a self-sustained oscillation 
(corresponding to a spiking state) is induced in the HH neuron model for J^c > /^c- Here, we 
set Idc = 100 jiK/cTo? and u to be the reciprocal of the golden mean [i.e., to = (v^ — l)/2], 
and numerically investigate dynamical responses of the (self-sustained) HH oscillator under 
the ac external stimulus. We first study the case of periodic forcing {i.e., A2 = 0) by 
varying Ai for /i = 26 Hz. Figures [T]^a)-[T|c) show the time series of V{t) for Ai = 50.42, 
50.33, and 50.24 /iA/cm^, respectively, and the bifurcation diagram in the Poincare map 
P is also given in Fig. [l|^d); stroboscopically sampled points in P are represented by solid 
circles in Figs. [T]^a)-[T|c). As Ai is decreased, successive period-doubling bifurcations occur. 
For example, periodic oscillations of V{t) in Figs. [l]^a) and[T|b) correspond to period-1 and 
period-2 states in P, respectively. When Ai passes a threshold A* (~ 50.28 /lA/cm^) a 
chaotic transition occurs. Thus, for Ai < A* chaotic oscillations with positive Lyapunov 
exponents appear, as shown in Fig. [l|c). 

From now on, we investigate the effect of quasiperiodic forcing on the period-doubling 
route to chaotic oscillation by changing Ai and A2 for /i = 26 Hz. Figure [2] shows a 
state diagram in the Ai — A^ plane. Each state is characterized by the largest (nontrivial) 
Lyapunov exponent a\, associated with dynamics of the variable x [besides the (trivial) 
zero exponent, related to the phase variable d of the quasiperiodic forcing] and the phase 
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sensitivity exponent S. The exponent S measures the sensitivity of the variable x with 
respect to the phase 6 of the quasiperiodic forcing and characterizes the strangeness of an 
attractor [12] . Regular quasiperiodic oscillations occur on smooth tori. A smooth torus that 
has a negative largest Lyapunov exponent {i.e., ai < 0) and has no phase sensitivity {i.e., 
6 = 0) exists in the region denoted by T and shown in light gray. When crossing a solid line, 
the smooth torus becomes unstable and bifurcates to a smooth doubled torus in the region 
represented by 2T. Smooth quadrupled tori, bifurcated from doubled tori, also exist in the 
region denoted by AT. On the other hand, chaotic oscillating states with positive largest 
Lyapunov exponents (cti > 0) exist in the region shown in black. Between these regular 
and chaotic regions, SN oscillating states that have negative largest Lyapunov exponents 
(cTi < 0) and positive phase sensitivity exponents {6 > 0) exist in the region shown in gray. 
Due to their high phase sensitivity, SN oscillating states have a strange fractal phase space 
structure. Various dynamical routes to SN oscillations via gradual fractalization, collision 
with a smooth unstable torus, and collision with a nonsmooth ring-shaped unstable set will 
be discussed below. 

When passing a heavy solid boundary curve in Fig. [2| a transition from a smooth torus 
to an SN attractor occurs via gradual fractalization [lOj. As an example, we study such 
transition to SN oscillations along the route a by decreasing Ai for = 0.1 /iA/cm^. Figures 
[3|^a)-|3|^c) show the time series of V{t) for the quasiperiodic oscillation, the SN oscillation, 
and the chaotic oscillation when Ai = 50.41, 50.374, and 50.36 /lA/cm^, respectively. These 
regular, SN, and chaotic states are analyzed in terms of the largest Lyapunov exponent ai 
and the phase sensitivity exponent 6 in the Poincare map. Projections of their corresponding 
attractors onto the 9 — V plane are shown in Figs. |3](d)-[3]^f). For the regular state, a smooth 
torus exists in the 9~V plane [see Fig. [3|^d)]. As Ai is decreased, the smooth torus becomes 
more and more wrinkled and transforms to an SN attractor without apparent mediation of 
any nearby unstable invariant set [3 [S]. As an example, see an SN attractor in Fig. ^e). 
This kind of gradual fractalization is the most common route to SN attractors. With further 
decrease in Ai, such an SN attractor turns into a chaotic attractor, as shown in Fig. Islf). 




A dynamical property of each attractor is characterized in terms of the largest Lyapunov 
exponent ai (measuring the degree of sensitivity to initial conditions). The Lyapunov- 



value of Ai ^ 50.377 /iA/cm^, an SN attractor appears. The graph of (Xi for the SN attractor 




exponent diagram ( 



i.e., plot of a I vs. Ai) is given in Fig. |3{g). When passing a threshold 
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is shown in black, and its value is negative as in the case of smooth torus. However, as Ai 
passes the chaotic transition point of Ai ^ 50.368 yuA/cm^, a chaotic attractor with a 
positive 0"! appears. Although SN and chaotic attractors are dynamically different, both of 
them have strange geometry. To characterize the strangeness of an attractor, we investigate 
the sensitivity of the attractor with respect to the phase 6 of the external quasiperiodic 
forcing [12]. This phase sensitivity may be characterized by differentiating x with respect to 
^ at a discrete time t = n. Using Eq. ([3]), we may obtain the following governing equation 
for ^ = 1,2,3,4), 

d fdxi\ dFi dxj dFi . , 

where (xi, X2, X3, X4) = {V,m,h,n) and Fj's (i = 1,2,3,4) are given in Eq. Starting 
from an initial point (x(0),^(0)) and an initial value d^/dO = for t = 0, we may obtain 
the derivative values of Sn'' (= dxi/dO) at all subsequent discrete time t = nhy integrating 
Eqs. (3) and (4). One can easily see the boundedness of Sn^ by looking only at the maximum 

7«(x(0), ^(0)) = max |5»(x(0), e(0))| = 1, 2, 3, 4). (5) 

Q<n<N 

We note that 7^^(x(0), ^(0)) depends on a particular trajectory. To obtain a "representative" 
quantity that is independent of a particular trajectory, we consider an ensemble of randomly 
chosen initial points {x(0),^(0)}, and take the minimum value of 7^'' with respect to the 
initial orbit points |12j . 

r^^ = . ±^nn^i^^WO)) (z = 1,2,3,4). (6) 

{x(0),6»{0)} 

Figure i^h) shows a phase sensitivity function rj^"*, which is obtained in an ensemble con- 
taining 20 random initial orbit points {(Vi(0), mj(0), /ij(0), ni(0), 6*1(0)); i = 1, . . . , 20} which 
are chosen with uniform probability in the range of Vi(0) G (—60,0), mj(0) G (0.1,0.9), 
hi{Q) G (0.1,0.2), miQ) G (0.5,0.7) and ^^(0) G [0,1). For the case of the smooth torus 
in Fig. i^d), F^^ grows up to the largest possible value of the derivative \dxi/d6\ along a 
trajectory and remains for all subsequent time. Thus, Fj^'* saturates for large and hence 
the smooth torus has no phase sensitivity {i.e., it has smooth geometry). On the other hand, 
for the case of the SN attractor in Fig. |3](e), T^^ grows unboundedly with the same power 
5, independently of 

f5;) ~ N^. (7) 



Here, the value of 5 ~ 2.39 is a quantitative characteristic of the phase sensitivity of the 
SN attractor, and 6 is called the phase sensitivity exponent. For obtaining satisfactory 
statistics, we consider 20 ensembles for each Ai, each of which contains 20 randomly chosen 
initial points and choose the average value of the 20 phase sensitivity exponents obtained 
in the 20 ensembles. Figure [3|^i) shows a plot of 5 versus AAi (= Ai — Al). Note that the 
value of S monotonically increases from zero as Ai is decreased away from the SN transition 
point Al (~ 50.377 /lA/cm^). As a result of this phase sensitivity, the SN oscillating state 
has strange fractal geometry leading to aperiodic complex spikings, as in the case of chaotic 
oscillations [e.g., see Figs. |3]^b) and[3]^c)]. 

As a dashed boundary curve in Fig. [2] is crossed, another route to SN attractors appears 
through collision between a stable smooth doubled torus and its unstable smooth parent 
torus [11]. As an example, we study this transition to SN oscillations along the route b by 
decreasing Ai for A2 = 0.06 /xA/cm^. Figure [4] shows a stable two-band torus (denoted by 
a solid curve) and an unstable smooth one-band parent torus (denoted by a short-dashed 
curve) for Ai = 50.348 /iA/cm^. The unstable parent torus is located in the middle of the 
two bands of the the stable torus. As Ai is decreased, the bands of the stable torus becomes 
more and more wrinkled, while the unstable torus remains smooth [see Fig. |4|^b)]. When Ai 
passes a threshold value of Ai ~ 50.3469 yuA/cm^, the two bands of the stable torus touch 
its unstable parent torus at a dense set of 9 values (not at all 6 values). As a result of this 
phase-dependent (nonsmooth) collision between the stable doubled torus and its unstable 
parent torus, an SN attractor is born, as shown in Fig. |4]^c). This SN attractor, containing 
the former bands of the torus as well as the unstable parent torus, has a positive phase 
sensitivity exponent {i.e., 6 > 0), inducing the strangeness of the SN attractor. However, 
its dynamics is nonchaotic because the largest Lyapunov exponent is negative {i.e., o\ < 0). 
As another threshold value of Ai ~ 50.344 /xA/cm^ is passed, the SN attractor transforms 
to a chaotic attractor with a positive largest Lyapunov exponent ai [see Fig. |4](d)]. 

A main interesting feature of the state diagram in Fig. [2] is the existence of "tongues" of 
quasiperiodic motion that penetrate into the chaotic region. The first-order (second-order) 
tongue lies near the terminal point (denoted by a cross) of the first-order (second-order) 
torus-doubling bifurcation curve. When crossing the upper boundary of the tongue (denoted 
by a dotted line), an intermittent SN attractor appears via phase-dependent collision of a 
stable torus with a nonsmooth ring-shaped unstable set [HI US] . We first study the transition 
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to an intermittent SN attractor along the route c in the first-order tongue by increasing A2 
for Ai = 50.34 /iA/cm^. Figure [S^a) shows a smooth torus for A2 = 0.093 /iA/cm^. When 
passing a threshold value of A2 — 0.093 35 yuA/cm^, a sudden transition to an intermittent SN 
attractor occurs, as shown in Fig. |5](b) for A2 = 0.093 53. Due to high phase sensitivity, this 
SN attractor with 5 ~ 3.15 has a strange fractal structure, while its dynamics is nonchaotic 
because of a negative largest Lyapunov exponent (ui ~ —0.015). A typical trajectory on 
the intermittent SN attractor spends a long stretch of time in the vicinity of the former 
torus, then it bursts out from this region and traces out a much larger fraction of the state 
space, and so on. To characterize the intermittent bursting, we use a small quantity d* for 
the threshold value of the magnitude of the deviation from the former torus. When the 
deviation is smaller (larger) than d*, the intermittent attractor is in the laminar (bursting) 
phase. For each A2, we follow a long trajectory until 10^ laminar phases are obtained in the 
Poincare map P and get the average of characteristic time r between bursts. As shown in 
Fig. |5]^c), the average value of r exhibits a power-law scaling behavior, 

r~AA^^, 7-0.5, (8) 

where the overbar represents time averaging and AA2 = A2 — A*2 {A*2 = 0.093 35). The 
scahng exponent 7 seems to be the same as that for the case of the quasiperiodically forced 
map pE]. As A2 passes another threshold value of A2 — 0.0936 /xA/cm^, the SN attractor 
transforms to a chaotic attractor because the largest Lyapunov exponent o"i becomes posi- 
tive, as shown in Fig. [sj^d). Furthermore, using the rational approximation, the mechanism 
for the intermittent route to SN attractors will be investigated below. Thus, a smooth torus 
is found to transform to an intermittent SN attractor via phase-dependent collision with a 
nonsmooth ring-shaped unstable set. 

We also study another intermittent route to SN attractors along the route d in the second- 
order tongue by increasing A2 for Ai = 50.3 /xA/cm^. Figure [sj^e) shows a smooth two-band 
torus for A2 = 0.03 /iA/cm^. As A2 passes a threshold value of A2 — 0.033 45 /lA/cm^, a 
band-merging transition from a smooth doubled torus to a single-band SN attractor occurs 
[T^ . Thus, an intermittent single-band SN attractor appears [e.g., see the intermittent SN 
attractor with ai ~ —0.029 and 5 — 2.17 in Fig. [5]^f) for A2 = 0.0336]. A typical trajectory 
of the second iterate of the Poincare map P {i.e. P^) spends a long stretch of time near one 
of the two former attractors {i.e., smooth tori), then it bursts out of this region and comes 
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close to the same or other former torus where it remains again for some time interval, and 
so on. As in the above case of intermittent route to SN attractors, we also obtain the 10^ 
laminar phases from a long trajectory in P^, and get the average of characteristic time r 
between bursts. As shown in Fig. [sj^g), the average characteristic time shows a power-law 
scaling behavior, 

Tr^AA'^, 7-0.5, (9) 

where AA2 = A2 — A2 {A2 = 0.033 45). The scaling exponent 7 seems to be the same as 
that for the case of the intermittent route to SN attractors occurring near the first-order 
tongue. Since the dynamical mechanism for the appearance of intermittent SN attractors 
near the first-order and second-order tongues are the same [i.e., an intermittent SN attractor 
appears via a phase-dependent collision between a smooth torus and a nonsmooth ring- 
shaped unstable set), the intermittent SN attractors for both cases seem to exhibit the 
same scaling behaviors. As A2 is further increased and passes another threshold value of 
A2 — 0.0352 yuA/cm^, the SN attractor turns into a chaotic attractor with a positive largest 
Lyapunov exponent ai, as shown in Fig. [Sj^h). 

The dynamical mechanisms for the appearance of intermittent SN attractors near the 
tongues are the same, irrespectively of the tongue order. Here, we consider the case of the 
main first-order tongue, and using the rational approximation to the quasiperiodic forcing, 
we search for an unstable orbit that causes the intermittent transition via collision with 
the smooth torus for A2 = 0.08 /iA/cm^. For the inverse golden mean co[= (-\/5 — l)/2], 
its rational approximants are given by the ratios of the Fibonacci numbers, Uk = Fk^i/F^, 
where the sequence of {F^} satisfies -Ffc+i = + -Ffc-i with Fq = and Fi = 1. Instead of 
the quasiperiodically forced system with u, periodically forced systems with ojk are studied 
in the rational approximation. As an example, we consider the rational approximation of 
level k = 7. The rational approximation to the smooth torus (denoted by a black curve), 
composed of stable orbits with period Fj(= 13), is shown in Fig. [6]^a) for Ai = 50.3432 
yuA/cm^. We note that a ring-shaped unstable set, composed of Fj small rings, lies near the 
smooth torus. At first, each ring is composed of the stable (shown in black) and unstable 
(shown in gray) orbits with period Fy [see the inset in Fig. [6|a)]. However, as Ai is changed 
such rings make evolution, as shown in Fig. ^h) for Ai = 50.343 /iA/cm^. For fixed values 
of Ai and A2, the phase 6 may be regarded as a "bifurcation parameter." As 6 is varied, a 
chaotic attractor appears via an infinite sequence of period-doubling bifurcations of stable 

10 



periodic orbits in each ring, and then it disappears via a boundary crisis when it colhdes 
with the unstable Fy-periodic orbit [see the inset in Fig. [6|^b)]. Thus, the attracting part 
(shown in black) of each ring is composed of the union of the originally stable F7-periodic 
attractor, the higher 2"F7-periodic {n = 1,2, . . .) and chaotic attractors born through the 
period-doubling cascade. On the other hand, the unstable part (shown in gray) of each ring 
consists of the union of the originally unstable Fy-periodic orbit [i.e., the lower gray line 
in the inset in Fig. ^h)] and the destabilized Fy-periodic orbit [i.e., the upper gray line in 
the inset in Fig. [6]^b)] via a period-doubling bifurcation. As the parameters, Ai and A2, are 
further changed, both the size and shape of the rings change, and eventually each ring is 
composed of a large unstable part (shown in gray) and a small attracting part (shown in 
black), as shown in Fig. |6](c) for Ai = 50.34 /lA/cm^ and A2 = 0.089 /iA/cm^ [a magnified 
view is given in Fig. [6|^d)]. We also note that new rings appear inside or outside the "old" 
rings. 

Finally, in terms of the rational approximation of level 7, we explain the mechanism 
for the intermittent transition occurring in the first-order tongue for Ai = 50.34 /iA/cm^ 
[see Figs. [5]^a)-[5]^b)]. As we approach the border of the intermittent transition in the state 
diagram, the ring-shaped unstable set comes closer to the smooth torus, as shown in Fig. |6](c) 
for A2 = 0.089 /iA/cm^ [see a magnified view in Fig. [6]^d)]. As A2 passes a threshold value 
of A2 — 0.090 305 /iA/cm^, a phase-dependent (nonsmooth) collision occurs between the 
smooth torus and the unstable part (shown in gray) of the nonsmooth ring-shaped unstable 
set. Then, the new attractor of the system contains the attracting part (shown in black) 
of the ring-shaped unstable set and becomes nonsmooth, which is shown in Fig. |6](e) for 
A2 = 0.0904 /iA/cm^ [see a magnified view in Fig. |6](f)]. As A2 is further increased, the 
chaotic component in the rational approximation to the attractor increases, and eventually 
for A2 — 0.091202 /xA/cm^, it becomes suddenly widened via an interior crisis when it 
collides with the nearest ring [e.g., see Fig.|6|^g)]. Then, Fj{= 13) "gaps," where no attractors 
with period F7 exist, are formed. A magnified gap is shown in Fig. [6|h). We note that 
this gap is filled by intermittent chaotic attractors. Thus, the rational approximation to 
the whole attractor consists of the union of the periodic and chaotic components. For this 
case, the periodic component is dominant, and hence the average largest Lyapunov exponent 
(((Ti) ^ —0.037) becomes negative, where (■ ■ ■ ) denotes the average over the whole 6. Hence, 
the rational approximation to the attractor becomes nonchaotic. We note that the 7th 
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rational approximation to the attractor in Fig. [6[g) resembles the (original) intermittent 
SN attractor in Fig. |5](b), although the level of the rational approximation is low. In this 
way, the intermittent transition to an SN attractor occurs through two steps in the rational 
approximation: the phase-dependent (nonsmooth) collision and the interior crisis. 



III. SUMMARY 



We have numerically studied dynamical responses of the quasiperiodically forced HH 
neural oscillator and compared them with those for the periodically forced case. For the 
case of periodic forcing, a transition from a periodic to a chaotic oscillation has been found 
to occur via period doublings in both numerical and experimental works. Effect of the 
quasiperiodic forcing on this period-doubling route to chaotic oscillation has been investi- 
gated. In contrast to the case of periodic forcing, new type of SN oscillating states have 
been found to exist between the regular and chaotic oscillating states as intermediate ones. 
Due to their strange geometry, these SN oscillations lead to the occurrence of aperiodic 
complex spikings, as in the case of chaotic oscillations. Hence, SN oscillating states might 
be a dynamical origin for the complex spikings which are usually observed in cortical neu- 
rons. Various routes to SN oscillations via fractalization, collision with a smooth unstable 
torus, and collision with a nonsmooth ring-shaped unstable set have been identified, as in 
the quasiperiodically forced logistic map [7] . These SN responses are also found to occur in 
other neurons exhibiting period-doubling route to chaos {e.g., the Morris-Lecar neuron and 
the FitzHugh-Nagumo neuron) under quasiperiodic stimulus [25]. Finally, we suggest an 
experiment on the quasiperiodically forced squid giant axon and expect that SN spikings to 
be observed. However, the real biological environment is a noisy one. Hence, it is necessary 
to further investigate the effect of noise on the SN response for real experiment. This type 
of investigation is beyond the scope of the the present paper, and hence it is left as a future 
work. 
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FIG. 1: Period-doubling route to chaos in the periodically forced HH oscillator for /^c = 100 
/iA/cm^ and /i = 26 Hz (A2 = 0). Time series of V{t) for (a) Ai = 50.42 ^A/cm^, (b) Ai = 
50.33 ^A/cm^, and (c) Ai = 50.24 fiA/cia^, which correspond to the period-1, period-2, and 
chaotic states in the Poincare map P (solid circles represent stroboscopically sampled points in P) , 
respectively. The largest Lyapunov exponent for the chaotic oscillation in (c) is ai c±. 0.096. (d) 
Bifurcation diagram {i.e., plot of V versus Ai) in P. 
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FIG. 2: State diagram in the Ai — A2 plane for i^c = 100 fiA/cia^ and /i = 26 Hz in the 
quasiperiodically forced HH oscillator. Regular, SN, and chaotic regions are shown in light gray, 
gray, and black, respectively. In the regular region, the torus, the doubled torus, and the quadrupled 
torus exist in the regions denoted by T, 2T, and AT, respectively, and the solid lines represent torus 
doubling bifurcation curves with terminal points denoted by crosses. SN attr actors appear via 
fractalization, collision with a smooth unstable torus, and collision with a nonsmooth ring-shaped 
unstable set when passing the heavy solid line, the dashed line, and the dotted line, respectively. 
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FIG. 3: Appearance of an SN attractor via fractalization along the route a in Fig. [2] for A2 = 
0.1 /xA/cm^. Time series of V{t) for (a) the quasiperiodic spiking state (^1 = 50.41 /iA/cm^), 
(b) the SN spiking state {Ai = 50.374 //A/cm^), and (c) the chaotic spiking state (Ai = 50.36 
/iA/cm^). Projections of attractors onto the 6 — V plane in the Poincare map are shown for (d) 
the smooth torus [corresponding to (a)], (e) the SN attractor [corresponding to (b)], and (f) the 
chaotic attractor [corresponding to (c)]. (g) Lyapunov-exponent diagram (i.e., plot of ai vs. Ai); 
a\ for the SN attractor is shown in black. (h)]^^hase sensitivity functions rj^^ are shown for the 
smooth torus (T) [(d)] and the SN attractor (SNA) [(e)]. For the case of the SN attractor, the 





FIG. 4: Appearance of an SN attractor via phase-dependent (nonsmooth) collision between a 
stable two-band torus and its unstable smooth parent torus along the route b in Fig. [2] for A2 = 
0.06 /iA/cm^. Stable two-band torus (denoted by a solid curve) and its unstable smooth torus 
(represented by a short-dashed curve) for (a) Ai = 50.348 /xA/cm^ and (b) Ai = 50.347 /iA/cm^. 
(c) SN attractor with ai ~ -0.015 and 6 ~ 3.77 for Ai = 50.346 fiA/cm^. (d) Chaotic attractor 
with cTi ~ 0.038 for Ai = 50.34 /xA/cm^. 
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FIG. 5: Appearance of intermittent SN attractors. (a)-(d) Appearance of an intermittent SN 
attractor along the route c in Fig. [2] for Ai = 50.34 /xA/cm^. (a) Smooth torus with ai ~ —0.071 
for A2 = 0.093 ixK/cvo?. (b) Intermittent SN attractor with o"i ~= —0.015 and 5 ~ 3.15 for 
A2 = 0.09353 /iA/cm^. (c) Plot of logioT vs. logioAyl2- The graph is well fitted with a dashed 
straight line with slope 7 ~ 0.5. Here r is the average characteristic time between bursts and 
A^2 = A2- AI {AI = 0.093 35 /iA/cm^). For each AA2, T is calculated from 10^ laminar phases 
in the Poincare map P. (d) Chaotic attractor with ai ~ 0.083 for A2 = 0.095 /^A/cm^. (e)-(h) 
Band-merging transition from a two-band torus to an intermittent single-band SN attractor along 
the route d in Fig. [2] for Ai = 50.3 ^A/cm^. (e) Smooth two-band torus with ai ~ —0.135 for 
A2 = 0.03 ^A/cm^. (f) Intermittent single-band SN attractor with cri ~ —0.029 and 5 ~ 2.17 for 
A2 = 0.0336 ;uA/cm^. (g) Plot of logior vs. logioAA2. The graph is well fitted with a dashed 
straight line with slope 7 ~ 0.5. Here r is the average characteristic time between bursts and 
A^2 = ^2 - ^2 (^2 = 0.03345 /iA/cm^). For each AA2, T is calculated from 10^ laminar phases 
in (h) Chaotic attractor with cji ~ 0.079 for A2 = 0.04 ^uA/cm^. 
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FIG. 6: Intermittent transition to an SN attractor in the rational approximation, (a)-(b) Birth 
and evolution of a nonsmooth ring-shaped unstable set for A2 = 0.08 /xA/cm^ in the rational 
approximation of level k = 7. (a) Smooth torus (denoted by a black curve) and ring-shaped 
unstable set for Ai = 50.3432 fiA/cm^. Each ring is composed of the attracting part (shown in 
black) and the unstable part (shown in gray). A magnified view of a ring is given in the inset, 
(b) Smooth torus (represented by a black curve) and an evolved ring-shaped unstable set for 
Ai = 50.343 /xA/cm^. As Ai is decreased, the unstable part (shown in gray) in each ring becomes 
dominant, (c)-(h) Mechanism for the intermittent transition to an SN attractor along the route c 
in Fig. [2]for Ai = 50.34 fiA/cra^. (c) Smooth torus and a ring-shaped unstable set for A2 = 0.089 
fiA/cra^ [a magnified view near 6 = 0.065 is given in (d)]. New rings appear inside or outside the 
old rings, (e) Nonsmooth attractor born via a phase-dependent (nonsmooth) collision between the 
smooth torus and a ring-shaped unstable set for A2 = 0.0904 ^uA/cm^ [(f) is a magnified view near 
9 = 0.065]. (g) Intermittent SN attractor with Fj gaps (born via an interior crisis) for A2 = 0.0915 
fiA/cTo^. (h) A magnified gap near 9 = 0.065. 
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